*************************************************************
* Replication codes for analyses in Section 3 of
* Beyond Yield Response: Weather Shocks and Crop Abandonment
* May 2020
*************************************************************	
		
	use dta\dta_for_reg, clear	

	rename acres_grain_ratio_harvest ha_corn
	rename acres_soy_ratio_harvest ha_soy
	rename l_yield_corn yld_corn
	rename l_yield_soy yld_soy
	
	***********************************************
	* Bin estimates 
	
	estimates clear	
	xtset fips year
			
	foreach x in corn soy {
		foreach y in ha yld {
	
			qui reghdfe `y'_`x' c.prec##c.prec v1_dd3bin* ///
				i.stateansi#c.t i.stateansi#c.t2 i.t if irr_ratio < .1 & drop_`x'==0 & pl_yrs_`x'>5, ///
				absorb(fips) vce(cluster fips stateyear)		
				estimates store `y'_`x'_dd3bin_2w	
			
			* obtain precipitation marginal effects (Figure 4)
				* margins, dydx(prec) at(prec = (1(0.2)16)) 
			
			qui xtreg `y'_`x' c.prec##c.prec v1_dd3bin* ///
				i.stateansi#c.t i.stateansi#c.t2 i.t if irr_ratio < .1 & drop_`x'==0 & pl_yrs_`x'>5, ///
				fe vce(cluster stateansi)		
				estimates store `y'_`x'_dd3bin_1w	
			
			}
		}
	
		* Results in Appendix Table A2 & Figure 2

			* two-way clustered results
			esttab *_dd3bin_2w, b(4) se(4) keep(*prec* *dd*) order(*dd* *prec*) staraux	
			
			* state-clustered results
			esttab *_dd3bin_1w, b(4) se(4) keep(*prec* *dd*) order(*dd* *prec*) bracket staraux		
			
		* Address LPM over-bound issue (Appendix Table A1)
			foreach x in corn soy {
				qui xtreg ha_`x' c.prec##c.prec v1_dd3bin* i.stateansi#c.t i.stateansi#c.t2 i.t if irr_ratio < .1 & drop_`x'==0 & pl_yrs_`x'>5, fe
				qui predict yhat_`x', xbu
				}
			sum yhat*, detail	
			
		* Bin robustness all 3C bins & all 4C bins (Appendix Figures A2 & A3)
		
		foreach x in corn soy {
			foreach y in ha yld {
				forvalues t = 1/3 {
					forvalues k = 3/4 {
		
					qui reghdfe `y'_`x' c.prec##c.prec v`t'_dd`k'bin* ///
						i.stateansi#c.t i.stateansi#c.t2 i.t if irr_ratio < .1 & drop_`x'==0 & pl_yrs_`x'>5, ///
						absorb(fips) vce(cluster fips stateyear)		
						estimates store `y'_`x'_v`t'_dd`k'bin
						}
					}
				}	
			}
			
		esttab *_dd3bin, b(4) se(4) keep(*prec* *dd*) order(*dd* *prec*)		
	
	***********************************************
	* Linear Spline Estimates (baseline and robustness)
		
		forvalues t = 25/27 {
			gen gdd8_`t' = (dday8C - dday`t'C)/100
			forvalues s = 34/36 {
				gen gdd`t'_`s' = dday`t'C - dday`s'C
				}
			}	
	
		foreach x in corn soy {
			forvalues t = 25/27 {
				forvalues s = 34/36 {
				qui reghdfe ha_`x' c.prec##c.prec gdd8_`t' gdd`t'_`s' dday`s'C ///
					i.stateansi#c.t i.stateansi#c.t2 i.t if irr_ratio < .1 & drop_`x'==0 & pl_yrs_`x'>5, ///
					absorb(fips) vce(cluster fips stateyear)		
					estimates store ha_`x'_dday`t'_`s'
					}
				}	
			}
			
		foreach x in corn soy {
			forvalues t = 25/27 {			
				qui reghdfe yld_`x' c.prec##c.prec gdd8_`t' dday`t'C ///
					i.stateansi#c.t i.stateansi#c.t2 i.t if irr_ratio < .1 & drop_`x'==0 & pl_yrs_`x'>5, ///
					absorb(fips) vce(cluster fips stateyear)		
					estimates store yld_`x'_dday`t'
					}
				}	
		
		* Present baseline estimates (Appendix Table A3)
		esttab ha*_dday26_35 yld*_dday26*, b(4) se(4) keep(*prec* *dd*) order(gdd8* gdd* dday* *prec*) staraux		
		esttab ha*_dday26_35 yld*_dday26*, b(4) se(4) keep(*prec* *dd*) order(gdd8* gdd* dday* *prec*)		
		
		* Effects evaluated at knots for ploting spline functions (Figure 3 & Appendix Figure A4)
		foreach x in corn soy {
			forvalues t = 25/27 {
				forvalues s = 34/36 {
					estimates restore ha_`x'_dday`t'_`s'
					lincom gdd8_`t'*(`t'-8)/100
					lincom gdd8_`t'*(`t'-8)/100 + gdd`t'_`s'*(`s'-`t')
					lincom gdd8_`t'*(`t'-8)/100 + gdd`t'_`s'*(`s'-`t') + dday`s'C*(41-`s')
					}
				}	
		}
		
		foreach x in corn soy {
			forvalues t = 25/27 {
				estimates restore yld_`x'_dday`t'
				lincom gdd8_`t'*(`t'-8)/100
				lincom gdd8_`t'*(`t'-8)/100 + dday`t'C*(41-`t')
				}	
		}		
			
		* state-cluster robust SE for baseline (Appendix Table A3)	
		foreach x in corn soy {
				qui xtreg ha_`x' c.prec##c.prec gdd8_26 gdd26_35 dday35 ///
					i.stateansi#c.t i.stateansi#c.t2 i.t if irr_ratio < .1 & drop_`x'==0 & pl_yrs_`x'>5, ///
					fe vce(cluster stateansi)		
					estimates store ha_`x'_dday26_1w
					
				qui xtreg yld_`x' c.prec##c.prec gdd8_26 dday26 ///
					i.stateansi#c.t i.stateansi#c.t2 i.t if irr_ratio < .1 & drop_`x'==0 & pl_yrs_`x'>5, ///
					fe vce(cluster stateansi)		
					estimates store yld_`x'_dday26_1w												
				}
		esttab ha*_dday*_1w yld*_dday*_1w, b(4) se(4) keep(*prec* *dd*) order(gdd8* gdd* dday* *prec*) bracket staraux		
						
	
	***************************
	* Two-season quad-avg spec
	
		foreach y in ha yld {
			foreach x in corn soy {	
				qui reghdfe `y'_`x' c.prec##c.prec c.tAvg_S1a##c.tAvg_S1a c.tAvg_S2a##c.tAvg_S2a ///
					i.stateansi#c.t i.stateansi#c.t2 i.t if irr_ratio < .1 & drop_`x'==0 & pl_yrs_`x'>5, ///
					absorb(fips) vce(cluster fips stateyear)		
					estimates store `y'_`x'_tAvg2
			}
		}	
		esttab ha*_tAvg2* yld*_tAvg2*, b(4) se(4) keep(*prec* *tAvg*) order(*tAvg* *prec*) staraux		
		
		foreach y in ha yld {
			foreach x in corn soy {	
				qui xtreg `y'_`x' c.prec##c.prec c.tAvg_S1a##c.tAvg_S1a c.tAvg_S2a##c.tAvg_S2a ///
					i.stateansi#c.t i.stateansi#c.t2 i.t if irr_ratio < .1 & drop_`x'==0 & pl_yrs_`x'>5, ///
					fe vce(cluster stateansi)		
					estimates store `y'_`x'_tAvg2_1w
			}
		}	
		
		* Results (Appendix Table A4)
		esttab ha*_tAvg2* yld*_tAvg2*, b(4) se(4) keep(*prec* *tAvg*) order(*tAvg* *prec*) staraux				
		esttab ha*_tAvg2_1w yld*_tAvg2_2w, b(4) se(4) keep(*prec* *tAvg*) order(*tAvg* *prec*) bracket staraux	
		
		* Marginal effects (Appendix Figure A5)
		estimates restore ha_corn_tAvg2
			margins, dydx(tAvg_S1a) at(tAvg_S1a = (0(4)32))
			margins, dydx(tAvg_S2a) at(tAvg_S2a = (0(4)32))
			
		estimates restore yld_corn_tAvg2
			margins, dydx(tAvg_S1a) at(tAvg_S1a = (0(4)32))
			margins, dydx(tAvg_S2a) at(tAvg_S2a = (0(4)32))					
		
	* One-season quad-avg spec

		foreach x in corn soy {
			foreach y in ha yld {	
				qui reghdfe `y'_`x' c.prec##c.prec c.tAvg##c.tAvg  ///
					i.stateansi#c.t i.stateansi#c.t2 i.t if irr_ratio < .1 & drop_`x'==0 & pl_yrs_`x'>5, ///
					absorb(fips) vce(cluster fips stateyear)		
					estimates store `y'_`x'_tAvg_1s
					}
				}
				
		esttab ha_corn_tAvg_1s*, b(4) se(4) keep(*prec* *tAvg*) order(*tAvg* *prec*) 
		
		* Marginal effects (Appendix Figure A6)
		
		estimates restore ha_corn_tAvg_1s
			margins, dydx(tAvg) at(tAvg = (0(4)32))
			
		estimates restore ha_soy_tAvg_1s
			margins, dydx(tAvg) at(tAvg = (0(4)32))			
			
		estimates restore yld_corn_tAvg_1s
			margins, dydx(tAvg) at(tAvg = (0(4)32))
			
		estimates restore yld_soy_tAvg_1s
			margins, dydx(tAvg) at(tAvg = (0(4)32))					

	
